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Abstract 

We propose a theoretical framework to study the cooperative behavior of dynamically coupled 
oscillators (DCOs) that possess dynamical interactions. Then, to understand synchronization phe- 
nomena in networks of interneurons which possess inhibitory interactions, we propose a DCO model 
with dynamics of interactions that tend to cause 180-degree phase lags. Employing an approach 
developed here, we demonstrate that although our model displays synchronization at high frequen- 
cies, it does not exhibit synchronization at low frequencies because this dynamical interaction does 
not cause a phase lag sufficiently large to cancel the effect of the inhibition. We interpret the disap- 
pearance of synchronization in our model with decreasing frequency as describing the breakdown 
of synchronization in the interneuron network of the CA1 area below the critical frequency of 20 
Hz. 

PACS numbers: 05.45.Xt, 75.10.Nr, 87.18.Sn 
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Many studies of cooperative behavior in many-body systems, of which one prototype 
is represented by spin systems, have been based on the assumption that the dynamics of 
interactions can be ignored. However, in the modeling of many biological systems, the 
inclusion of interaction dynamics is necessary to obtain a faithful description. In much of 
the research that has been done, systems with such slow interaction have been mimicked by 
a so-called multi- compartment model which corresponds to a discrete version of a distributed 
parameter system. By using such a multi-compartment model, we can accurately estimate 
the dynamical behavior of a single unit. However, it is impossible to carry out a large-scale 
numerical simulation with such models to understand the behavior of a large-scale system 
composed of tens of thousands of units. Even if it is possible, we cannot understand the 
essence of such phenomena. To address this problem, we need to reduce the degrees of 
freedom needed in a descriptive model. 

Phase reduction is a powerful tool that helps us understand the oscillatory behavior 
of large-scale systems composed of such multi-compartment models. In previous studies 
of systems consisting of such multi-compartment models, phase descriptions of the entire 
systems including the dynamical interactions themselves were obtained numerically [[I], |2|, |3|, 
Q. Thus, the effects of the interaction dynamics were numerically incorporated into only a 
phase- coupling function to establish the phase equation. Therefore, in general, the effects of 
the interaction dynamics alone cannot be extracted in such methods. 

In this paper, we analyze systems of dynamically coupled oscillators (DCOs), and demon- 
strate how the dynamics of the interactions influence the cooperative behavior exhibited by 
the systems. Using the multi-scale perturbation method (MSPM) || |6| to carry out a phase 
reduction, we are able to represent the nonlinear dynamical interaction in terms of only the 
amplitude and phase of the fundamental frequency component in the frequency response 
of the interaction. In control engineering, the amplitude and phase of the fundamental 
frequency component of the frequency response are referred to as the describing function. 
The describing function used in the treatment of non-linear systems is an extension of the 
transfer function used in the treatment of linear systems. In contrast to previous work 

[2|, H |J, the effects of the interaction dynamics are expressed by the describing function 
in our approach, and so it is easy to isolate and extract the effects of the interaction dy- 
namics from the phase equation. Therefore, our approach can clearly elucidate the effect of 
the interaction dynamics on the cooperative behavior of entire systems. Our results imply 
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the possibility of a new modeling approach for biological systems by which the describing 
function of the interaction is identified. The description proposed here corresponds more to 
higher order approximation than to conventional linear description (e.g., the auto-regressive 
(AR) model). 

In the first half of this paper, we show in detail the process of analytical phase reduction 
for a general model with dynamical interaction. We treat a large population of Stuart- 
Landau (SL) oscillators coupled through nonlinear dynamical interaction. The SL oscillators 
have the essential structure of the Hopf-bifurcation, because its evolution equation can be 
derived from any non-linear oscillator system with the Hopf-bifurcation through perturbation 
expansion ||. In the latter half of the paper, we apply our theoretical framework to the 
analysis of a network of interneurons. The analysis of an interneuron network is a good 
example to verify the applicability of our approach for understanding the neural cooperative 
behavior via dynamical interaction. 

First, we consider a large population of SL oscillators weakly coupled through interac- 
tions which themselves are nonlinear dynamical systems. We analyze the following coupled 
system: 



Here Wj is the state variable (a complex number) of the jth oscillator (with a total of 
N). Q is the average natural frequency and eUj represents the deviation from the average 
for the jth oscillator randomly distributed with a density represented by g(u). The terms 
multiplied by e are considered perturbations, and thus, the quantity e controls the magnitude 
of perturbations. When e = 0, the system has an unstable fixed point at the origin and a 
stable limit-cycle orbit on the unit circle in the complex plane represented by 



where <pj is the phase of the jth oscillator. When e = 0, the system is neutrally stable with 
respect to a perturbation in the form of a temporal shift while it conserves a fixed orbit. 




(1) 



k=l 



wj(t) = $j(t), $j(t) = expi (fit + (j>j) , v 0j G R/2tt, 



(2) 



Thus, <f>j depends only on the initial condition. The quantity Sjk in Eq. ([I]) represents the 
output from the dynamical system of the interaction expressed by 




(3) 
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where Xjk represents the internal state of the interaction, Fjf-(-) is a nonlinear function 
determining its dynamics, and Vjk(-) is an output function. 

We now discuss the dynamical interaction described by Eq. (|3]). If Wk = = 
exp i (Qt + 4>k) is input into Eq. (|3p, the higher harmonics resulting from the nonlinearity of 
this equation are superimposed on Sjk- Here we restrict our consideration to the case that 
the output Sjk is a periodic function possessing the same period as the input Wk- In this 
case, Sjk can be expanded into the following Fourier form: Sjk = jfk^k + Yl™>i Jjk ^ 
Eq. (H) takes a form of a linear system, jjj^ would be non-zero, and jj^ for n > 1 would 
all be zero. In this case, J-^ is called the transfer function of the dynamical interaction Eq. 
(0). If, on the other hand, Eq. takes a form of a nonlinear system, jj^ and some of J-^ 
for n > 1 would be non-zero. In this case, jj], is called the describing function in the field 
of control engineering. 

Now, employing the MSPM, we attempt to reduce Eq. (P to a phase equation. The 
solution of the perturbed system ([]]) can be represented as 

Wj (t) = $j(t) + euj(t), $j(t) = exp i (fit + <j>j(et)) , (4) 

where we write <pj as a function of et to make explicit its slow evolution in time, and euj 
is the first order change in Wj introduced by the perturbation. If we substitute Eq. (|J) 
into Eqs. (|l|) and (|j) and expand the resulting expression about e = 0, the 0(e) equation 
becomes 

^^t) = ^ + L ^ + if £ £ 4 n) ^- ( 5 ) 

^ ' n=l fc=l 

Here is the linear operator corresponding to Eq. ([!]) linearized about the periodic solution 
(@) in the case e = 0. It is given by 

L^u = + ^fiw ~ 3* 2 ^ ~ M ; (6) 

where the overline denotes complex conjugation. All of the eigenvalues of are non-positive, 
since the solution $(£) is stable. Note that there exists the eigenfunction $'(£) of with 
eigenvalue 0. This eigenfunction corresponds to an infinitesimal temporal shift because 
<3?(t + 5) ~ $(£) + £$'(£). We assume there exists no other eigenfunction with eigenvalue 
in the space of the periodic functions, so that kerL^, = span{$'(t)} = span{i$(t)}. This 
assumption is equivalent to that of the orbital stability of 
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We define the inner product of two ^y-periodic complex functions u(t) and v(t) as 
(u(t),v(t)) = jfi dt(Re{u(t)}Re{v(t)} + lm{u(t)}lm{v(t)}). Then, we obtain the oper- 
ator adjoint to as 

Liu = ^ - itlu - $ 2 u - u. (7) 
v dt w 

From Fredholm's alternative 0, there is a so-called response function $*(£) that spans the 
kernel of L£ in the space of periodic functions. Therefore, kerL^ = span{$*(t)}. In this 
case, we explicitly obtain = Taking the inner product of i$(t) and Eq. (|^), we 

obtain 

(< ** ^ ~ ^ u >-^ 

k,n 

= L^Uj) = (/,',/'I , ,.<' ; ) = 0, (8) 

where the quantity et is treated as a constant over a single period, since <f)j{et) is driven by 
a weak perturbation. Here, the higher harmonics resulting from the nonlinearity of Eq. (|3|) 
cancel, because jfk^l^j = (n > 1). Thus, we derive the phase equation describing 

the slow phase dynamics of the system, 



^Jy = Uj + A i k sin ^ k - fa + M > ( 9 ) 

j 

with A jk = \J$\, ip jk = arg jf k \ (10) 



where Uj represents a natural frequency of unit j at the level of the phase equation. As 
Eqs. (H) and ([H]) reveal, in the phase equation, the dynamical interaction (|3|) is expressed 
in terms of the describing function. Thus, seen the fundamental frequency components in 
the output from the dynamical system of the interaction (^) clearly play a key role in the 
slow phase dynamics of the system. 

If all interactions are identical, i.e. J-^ = J (Aj k = A and ijjj k = ip), Eq. @ is equivalent 
to the Sakaguchi-Kuramoto model [0], and so, in this case, the Sakaguchi-Kuramoto theory 
can be applied to the analysis of Eq. (H). Here, we define the mean field asm = | Sj=i w j- 
When \m\ ^ 0, the system is in a state of phase synchronization. In the thermodynamic 
limit, N — > oo, we obtain the following self-consistent equation relating the order parameters 
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FIG. 1: (a) Dynamical behavior of neural components, (b) Schematic diagram of the dynamical 
interaction employed in our model which consists of a rectangular-pulse generator and a second- 
order lag system (SOLS), (c) Bode diagram of the SOLS, G(ifl) = 1/(1 + iOr) 2 . 
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Here, to derive the order parameter equation, we have assumed the existence of one large 
cluster of oscillators phase-locked at frequency Cl in the slow dynamics described by ([5]). In 
the original system, described by flip, the frequency of the synchronous cluster is Q + efl. 

Next, we apply our theoretical framework to the analysis of a network of interneurons as 
an example. High frequency oscillations, known as 7 oscillations, have been observed in many 
areas of the brain. The synchronization of inhibitory interneurons is known to be the origin 
of 7 oscillations ||. Inhibitory interactions, which correspond to anti-ferromagnetic inter- 
actions in spin systems, introduce competition; for this reason, in general, such interactions 
tend to prevent synchronization. Therefore, the synchronization phenomena in networks of 
interneurons are non-trivial [Q, || |10| . A series of physiological in vitro experiments inspired 
our study. One in particular was the experimental discovery that the interneuronal network 



in the CAl area of the hippocampus displays synchronization at frequencies greater than or 
close to 20 Hz, but not at frequencies significantly below this value ]TT|, [12[ . The discovery 
of this phenomenon provides evidence suggesting that cooperative behavior in networks of 
interneurons arises due to the dynamical nature of the interactions. 

To understand the mechanism of the synchronization, we propose a DCO model with 
interaction dynamics that tend to cause 180-degree phase lags. The estimated connectivity 
between interneurons in the CAl area is 10% UTO"] . The interneuron network in the CAl 
area has dense connectivity, because the number of connections between interneurons is 
O(N). Therefore, we can capture properties of the interneuron network by analyzing the 
global coupling systems. As shown in Fig. |l](b), the dynamical interactions of neurons are 
modeled by a rectangular-pulse generator consisting of a threshold element and a second- 
order lag system (SOLS): 
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The parameters rj^ and are the time constants of the two cascade units of the first-order 
lag system that form the SOLS. In Eq. fljjD, the parameters G and H denote, respectively, 
the gain and the threshold of the threshold element, and are defined as G = 1/6 and 
H = cos 6 (6 > 0). If the value of 9 is small, the rectangular pulse generated by the 
threshold element is sharp. Note that Eq. fll2"|) represents mutual interactions only between 
the real parts of the SL oscillators. 

In the context of control engineering, the SOLS is regarded as one of the minimal models 
for realizing a 180-degree phase lag. As shown in Fig. 0(a), the electrical properties of somata 
and proximal dendrites can be approximated as the CR circuit. Near resting potential, we 
can obtain a linear approximation of their behavior as being that of a first-order lag system. 
The response of the synaptic conductance is also displayed in Fig. [l|(a). The evolution of 
the synaptic conductance is roughly reproduced by that of the first-order lag system. Thus, 
as shown in Fig. [l|(a), the SOLS given by Eq. fljjD can be considered as a cascade system 
consisting of these components. In addition, to model the generation of action potentials in a 
neuron when its membrane potential reaches the threshold, and furthermore, to demonstrate 
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FIG. 2: (a) \m\ as a function of a. The data points with error bars represent the means and standard 
deviations of 10 samples obtained from numerical simulations of Eq. ([l]). Here, $7 = 2-7r80, t = 0.02, 
9 = vr/6, e = 10 and N = 1000. (b) Critical deviation a c as a function of fir. The quantity plotted 
here is actually a normalized critical deviation, with maximum value 1. 

the applicability of our approach in analyses of non-linear systems, we have also introduced 
a threshold element into the dynamical interaction given in Eq. (|I2"D. 
The describing function for the dynamical interaction given in flT2"|) is 

2sin# 1 



J 



(i) 



« 9 (l + <r«n)(l+ir^n) 

Here, the describing function of the threshold element is 2 sin#/ (tt9), which contributes only 
to the gain of Eq. ([13]). Figure |l|(c) displays the Bode diagrams of the SOLS. As Fig. |l](c) 
reveals, the gain of the SOLS, which is a type of low-pass filter, is a decreasing function 
of frequency at high frequencies, while the phase of the SOLS converges to -180 degrees as 
the frequency increases. If the phase lag is 180 degrees, J-j^ is a positive real number, and 
the system thus effectively becomes a ferromagnetic oscillator system. Note that Eq. (|I2"D 
represents the mutual interactions only between the real parts of the SL oscillators. Thus, 
Ajk takes half the value it had in Eq. fllPf) , and here Ajk = §1^1, ipjk = arg J.^ . 

In the numerical calculation we will now discuss, we chose the distribution of natural 
frequencies as g(u) = (27ro" 2 )~ 1//2 exp (— u 2 /(2a 2 )), where a is the deviation of the natural 



_(2), 



(13) 



frequencies. We begin with systems where the interactions are all equivalent, i.e., fir- 



(i) 



Q, = fir, and hence Ajk = A and ipjk = if). In this case, Eq. @ is equivalent to the 
Sakaguchi-Kuramoto model 0. Figure 0(a) displays \m\ as a function of o in the case of 
the dynamical interaction given by Eq. (0), where the solid curves were obtained from 
the order parameter equation (|TT|), and the data points with error bars represent results 
obtained using the fourth-order Runge-Kutta method with Eqs. ([!]) and (0). Next, we 
estimate the critical value of a, representing the boundary between \m\ = and \m\ ^ 0. 
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FIG. 3: System size scaling. Deviation of local fields multiplied by v N as a function of N. 

Figure |2](b) displays the critical deviation a c as a function of Qt for the two models. Note 
that the absolute values of A and a c have no meaning in comparing such models, because 
a c is proportional to A 0. Plotted in Fig. ^|(b) is the critical deviation normalized to have 
a maximum value of 1. The widths of the synchronous regions decrease as the frequency 
increases because, as shown in Fig. 0(c), the SOLS is a low-pass filter whose gain rapidly 
decreases at higher frequencies. However, if G in Eq. ([12]) is sufficiently large, the gain 
of the interaction given by Eq. ( |L2|) can be made sufficiently large by the pulse generator 
to realize synchronization at high frequencies. However, as Fig. ||](b) reveals, our model 
does not exhibit synchronization at low frequencies because, as shown in Fig. [l|(c), this 
dynamical interaction does not cause a phase lag sufficiently large to cancel the effect of the 
inhibition. We interpret the disappearance of synchronization in our model with decreasing 
frequency as describing the breakdown of synchronization in the interneuron network of the 



CAl area below the critical frequency of 20 Hz pi] , |12| . Based on this correspondence, we 
conjecture that in an interneuron network, a phase lag resulting from the dynamical nature 
of the interactions might cancel the effect of inhibition at high frequencies, and through this 
cancellation, the system effectively becomes a ferromagnetic oscillator system. 

By means of numerical simulations, as shown in Fig. [|, we have demonstrated through 
a scaling plot that the variation of a local field is 0(1/ \/N) in systems where a quenched 
random time constant is used for the interaction dynamics in Eq. (|12"1). In the thermodynam- 



ical limit iV — > oo, the system asymptotically approaches a system with global homogeneous 
interaction, even with a quenched random time constant of the interaction dynamics. There- 
fore, as the system size N increases, the system tends to synchronize through the cancellation 
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of fluctuation via interaction. 
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